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ABSTRACT 

Aims. General 3D astrophysical atmospheres will have random velocity fields. We seek to combine the methods we have developed 
for solving the ID problem with arbitrary flows to those that we have developed for solving the fully 3D relativistic radiative transfer 
problem for monotonic flows. 

Methods. The methods developed in the case of 3D atmospheres with monotonic flows, solving the fully relativistic problem 
along curves defined by an affine parameter, are very flexible and can be extended to the case of arbitrary velocity fields in 3D. 
Simultaneously, the techniques we developed for treating the ID problem with arbitrary velocity fields are easily adapted to the 3D 
problem. 

Results. The algorithm we present can be used to solve 3D radiative transfer problems that include arbitrary wavelength couplings. 
We use a quasi-analytic formal solution of the radiative transfer equation that significantly improves the overall computation speed. 
We show that the approximate lambda operator developed in previous work gives good convergence, even neglecting wavelength 
coupling. Ng acceleration also gives good results. We present tests that are of similar resolution to what has been presented using 
Monte-Carlo techniques, thus our methods will be applicable to problems outside of our test setup. Additional domain decomposition 
parallelization strategies will be explored in future work. 
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1. Introduction 



The photon's four-momentum can be written 



In a series of pa pers ([Hauschildt & Baron] |2006| jBaron & 
"Hauschildt 2007; Hauschildt & Bai'on 2008 2009 ; Baron et al. 



|2009b,, Hauschildt & Baron,2010,.Seelmann et al. 2010,) we have 
developed a general framework for solving 3D radiative trans- 
fer problems in Cartesian, cylindrical, and spherical coordinates 
for both static and monotonic velocity fields in the comoving 
frame. We have also developed a n Eulerian code for velocities 
< 1000 km s ' (Seelmann et al. 2010). The neglect of relativis- 
tic effects and resolution constramts limits the applicability of 
the Eulerian approach to v/c < 0.01. Our affine method for solv- 
ing the fully relativistic transfer equation is exact to all orders 
in v/c ( |Chen et al.||2007| [Baron et al.1|2009a| l. In a parallel se- 
ries of papers we have introduced methods for solving the fully 
relativistic transfer equation in ID spherical coordinates with ar- 
bitrary v elocity fields ( [Baron & Hauschildt|[2004[ [Knop et al] 
|2009a|bD . 



2. Basic formalism 



Using the general formalism developed in lChen et al.[ ( [20 07^ we 
can derive the transfer equation in flat spacetime with arbitrary 
flows. We choose to work in spherical coordinates without loss 
of generality. 



a dx" h 

/loo 



(1) 



where h is Planck's constant, ^ is the affine parameter, /loo is the 
rest frame wavelength, and n is the 3D direction of the photon as 
seen by a distant, stationary, observer. The four-velocity of the 
comoving observer in an arbitrary flow can be written 

=r(r,0[l,)8(r,f)], (2) 
and the comoving wavelength A can be obtained using 

7 = -(« ■ P)- (3) 

A 

Here, jS = v/c, and j - {\ are the usual quantities of 

special relativity. The 3D geodesic in the flat spacetime can be 
parametrized as 

r(i) = ro + ft (4) 

where ro is the starting point of the characteristics, and s is the 
rest frame physical distance related to the affine parameter ^ by 



h 

/too 



(5) 



1 



E. Baron et al.: Arbitrary Velocity Flows 



The radiative transfer equation can be written in terms of the 
affine parameter ^ as (see Eq. (10) of Chen et al.|2007 1: 



When we numerically integrate the radiation transfer equation, 
a(s) can be approximated as 



dA dlx 



h 5dA 
^^A ^ Ad^ 



h 



(6) 



a{s) 



6s{\-ri-l3) 



(16) 



where /^i(r, f; n) is the specific intensity measured by a comoving 
observer (note that 7^/1^ is observer-independent) at the (global) 
rest frame space-time point (r, f), toward the rest frame direction 
n, and at the comoving wavelength A. When expressing the 7-D 
phase-space dependence of the comoving specific intensity, the 
only comoving variable we used was the wavelength A, in par- 
ticular, we did not use angles measured by a comoving observer 
to specify the direction of the photons. The advantages for this 
at first sight odd phase-space configuration have been explained 
in detail in Chen et a"Ll ( |2007 1, and we applied this technique in 
[Baron et al.| ( |2009a| l. 

We can rewrite Eq. (|6]l as 



where 6s is the differential step size (physical distance) along 
the characteristics, 6{n ■ P) and Sfi are the changes ofh-fi and /3, 
respectively, when we move one step forward which includes the 
changes induced by both time and spatial advances, for instance 



(17) 



Since few numerical schemes will be able to provide the fully 
implicit derivative, 6/3, will often be obtained for example by the 
backward difference 



6/3 = {/3(si,ti)-/3(si,ti^i) + (J3(si,ti) -/3(si^utd- 



(18) 



d{ct) 1 dh 
~dfc~dt 



+ — V/^ 

d^ 



dAdh 
d^'dA 



In the stationary case, both (3 and f(s) are independent of 
time and specializing Eqs. Q and (12 1 to that case, a(s) becomes 



h 5dA\^ h 
-^'^'A^Ad-^r^'^'A' 



a(s) — 



n ■ V(n ■ ^S) - y^/3(l - n ■ j8)(n ■ V/3) 



or equivalently 



1 -nj8 

= -(n ■ V) ln(l - n ■ yS) - y^/3(n ■ V/3), 



(19) 



d{ct) 1 diA 
~dfc~dF 



ds dr ds dA dl,i 

d^ ds d^ ds dA 



-\Xa 



h 5 ds dA\ h 

^Ad-^dsy^^'^^^^ 



A A d^ ds I A 
Then using the definition of s from Eq. (jSj) and the fact that 



where we have used the fact that along the characteristics, d/ds 
no longer contains the time derivative and is thus the directional 
derivative operator, that is, d/ds = fi ■ V. Recall that in the 
flat spacetime that we are considering, our characteristics are 
straight lines for all velocity flows. 

In terms of its spherical components, /3 can be written 



P^/3rtr +/3eh +/3^t^ 



(20) 



d{ct) 
~d^ 



Ao, 



(9) 



where er,efl,e0 are the spherical orthonormal basis vectors at 
point r(r, 0, </>), i.e. 



we find 

1 dl,i 
elk 



A ds 



dA dlx 
i ds dA 



/Im 5 dA\ 

= -\^'^^-aTsV' 



Cr = (sinflcos^, sin0sin0, COS0), 
Cfl - (cos 0COS 0, cos 0sin(/>, - sin0), 
- (- sin (p, cos (p, 0), 



(21) 



(10) 



and consequently the n - P 'm Eq. ( 16 1 can be calculated using 



Eq. ( 10 1 can be put into our standard form: 



/3rn.-tr+/3ga-te+) 
= /3rnr +/3eng + /3^n^ 



(22) 



Idh 
c dt 

where 



dl d 
^ + a{s)—{Ah) + Aa{s)h = -XAf{s)lA + riJis), 



ds 



dA 



r(r,0[l-fl-j8(r,f)] 



is simply the Doppler factor, and 



a{s) = 



1 dA 
A ds 



Using Eqs. Q and ( [T2] i, a(s) is found to be 
1 



a(s) 



1 



d . 

d-s^^-^' 



n/3 

where /? is the magnitude of fi, and 



7^/3(1 



h/3) 



dl 
ds 



(11) 



(12) 



(13) 



(14) 



(note that along the characteristics, n has constant Cartesian 
components, n^, n,, but changing spherical components, 
rir, rig, n^). Writing n = (n^-,ny,n-), the spherical components 



ng, «0 can be easily computed using Eq. (21 



2.1. Comparison with Mihalas 



At first glance comparing Eq (111 with Eq. (2.12) of Mihalas 
( 1980 1 something seems amiss. Like Mihalas ( |1980 1, we work in 
the frame where spatial coordinates and clocks are measured by 
an observer at rest. However, Mihalas' time derivative contains 
a Doppler factor, whereas ours does not. Also, our terms on the 
right-hand side contain Doppler factors, f(s), whereas those of 
Mihalas do not. The discrepency has been noted in passing by 
Che n et al.| (j2007) and arises because our 5 is a true distance 
measured in the observer's frame, whereas that of Mihalas, sm, 
contains an extra Doppler factor: 
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(note that fis) is not a constant along the characteristics, and 
therefore is not related to the physical distance s by a simple 
afiine transformation). Thus, we can transform from s to sm in 
Eq. ([T0| to find 



A 1 dh 

/loo c dt 



+ ^ 



dA dh 
X dsM dA 



I 5 dA \ ^ ^^^^ 
\XA + ^:i—\h + r]A, (23) 



which is very similar to the equation of Mihalas, except that 
the coefficient multiplying the time derivative term is the inverse 
Doppler factor /(i)"' , since we are working with instead of /,„ 



as did Mihalas. Jack et al. (2009 1 accidentally forgot to convert 
the time derivative from ly to and thus their Eqns. (24)-(25) 
have a coefficient of the time derivative with the Doppler factor 
in the numerator, rather than in the denominator. Thus, the time 
derivative terms in their Eqns. (24)-(25) should be multiplied by 
fis)"^ to give the correct equations. 



2.2. Nonmonotonic flows 



We are now in a position to tie together the work of Baron 
& Hauschildt] (p004|, |Knop et al.] (|2009a) and 'Baron et al 



( |2009b[ ). The formal solution in the monotonic case is an ini 
tial value problem in wavelength, but in the arbitrary flow case 
both spatial coordinates and wavelengths m'e fully coupled. This 
poses a significant memory cost, since the matrix obtained by 
finite-diffferencing the equations now contains every wavelength 
and not just the spatial points along the characteristic. The 
computational cost is surprisingly low because the linear sys- 
tem can be solved using the semi-analytic method of Knop 
et al. ( 2009a[l. Wh ile the framework that was given in Baron & 



Hauschildt| (|2004) andl Baron et al.| ( |2009"F ) was formulated for 



just one wavelength discretization, we included here the fully 
implicit discretization developed in Hauschildt & Baron ( 2004) l. 
Furthermore, we used the new formal solution that avoids nega- 
tive generalized opacities ( |Knop et al.|2009b} . 

The stationary equation of radiative transfer in its character- 
istic form for the specific intensity / along a path s reads 



^ = m^li - f{s)xih - 4a,I, - a, ^ , (24) 
as oA 



where 77 is the emissivity, x the opacity, and the subscript I in- 
dicates dependence on wavelength. The J^-term is the coupling 
term between the wavelengths and depends on the structure of 
the atmosphere and on the mechanism of the coupling (jMihalas 
[T980l l. 

The wavelength derivative can be discretized in two ways 
as described in [Hauschildt & Baron| (2004). The different dis- 
cretizations can be mixed via a Crank-Nicholson-like scheme 
with a mixing parameter ^ e [0, 1]. The wavelength discretized 
equation of radiative transfer can then be written as 



As 



f{s)m-f{s)xiIi-a,[A + ^l))li 

-^a,{pjli-x + ptii+i) 

-[1 - ^] a, [pjli-i + p\li + pUm) , 



where the p' coefficients in an ordered wavelength grid < 
Al < Ai+i are defined as 



Pi 
P\ 

p! 
Pl 
p\ 

pI 



Al- 



1 

Al 


- Ai-i 




Al 


^1 


-A,-i 














Al 


Al 


- Ai+i 




Ai+\ 


Al 


- Ai^\ 



for at, > 



for ua, < 0. 



(26) 



(27) 



The dependence on the sign of oa is introduced to define local 
upwind schemes (see Baron & Hauschildt|2004| l. 



After introducing a generalized opacity (see Knop et al. 
2009a.b) 

XI ^ mxi + ^ aip\, (28) 
defining the source functions 



Si = 
Si = 



(29) 
(30) 



fis)Si-^^^{pJIl-l+pllM)'^ 

Si = (/'r'f/-i+K//+i) + (4 + [l-^];^l)/j, (31) 

a formal solution of the radiative transfer problem can be formu- 
lated. We used a full characteiistic method throughout the atmo- 
sphere. The spatial position on a characteristic is then discretized 
on a spatial grid. In the following a pair of subscript indices will 
mark the position in the spatial grid and in the wavelength grid. 
Commonly the spatial grid is mapped locally onto an optical 
depth grid via the relation dr/ - xids. The formal solution of 
the radiative transfer equation (25 1 between two points and 
Sj on a spatial grid along the photon path can be written in terms 
of the optical depth as follows: 



6hi 



= \ Sie^ ^'dr = a,jS + /3ijS + yijS 
]j = I Sie''''''dT^aijSi-ij+$ijSij, 



i+\,i 



(32) 
(33) 

(34) 



with At - T,+i ; -T/ / and T/ / - f '' ;^;(^)d^. The a-jS-y coefficients 
are described in [Olson & Kunasz| ( fTM7l ) and [HauschIIdt| ( [ 1 992| l . 
dii in Eq. (34 1 is linearly interpolated and in general difi'erent 
from the coefficients in Eq. ( [33] l and is therefore marked with a 
tilde. 



(25) 



Eq. ( 32 1 can be wiitten in matrix notation for any given char- 
acteristic: 

I = A ■ I + Al. (35) 

Here I is a vector with all intensities, A is a square matrix that de- 
scribes the influence of the different intensities upon each other, 
and Al is a vector with the thermal emission and scattering con- 
tribution of the source function. For a characteristic with n, spa- 
tial points and «/ points in the wavelength grid, the intensity vec- 
tor I has n, X n; entries. In the following a superscript of k will 
label the characteristic at hand. The components of the matrix A 
from Eq. ( 35 1 at the spatial point / and the wavelength point I are 
given by 
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of the radiative transfer equation, which is symbolically written 
using the A-operator as 



. k "/+!,; -,k 



X, 



Xi-ij 



B 



c: 



a' 



B1 



+,k 



Xi-ij 



C 



+,k 



k 1+1,/ +,/t 



(36) 

(37) 

(38) 

(39) 

(40) 
(41) 
(42) 

(43) 

(44) 



J A - a- 

For the transition of a two-level atom, we have 



(47) 



(48) 



where J - j ^{^)Ja dA, A = J (/)(A)A^ dA with the normalized 
line profile (j)(A). 



The A-iteration method, i.e. solving Eq. (48 i by a fixed-point 
iteration scheme of the form 



new — AS r 



'5 new — (1 f)-^new + ^B, 



(49) 



fails in the case of high optical depths and small e. This is 
because the largest eigenvalue of the amplification matrix (for 
Doppler-profiles) is approximately /Imax ~ (l-e)(l-r~'), where 
T is the optical thickness of the medium ( [Mihalas et al.|1975] l. 
For small e and high T, this is very close to unity and, therefore, 
the convergence rate of the A-iteration is very poor. A physical 
description of this effect can be found in Mihalas | ( |1 9 80| l. 

The idea of the ALI or operator splitting method is to re- 
duce the eig envalues of the amplification matrix in the itera- 
tion scheme ( |Cannon||l973 1 by introducing an approximate A- 
operator (ALO) A* and to split A according to 



A = A* -H (A - A*) 



Following Knop et al.l(l2009al), the naming scheme of the quan- 
tities in Eqs. (36l-(44i indicates with which specific intensity 
element they are associated. For an index pair / and / a super- 
script refers to an intensity at wavelength I - 1, a superscript 
to the same wavelength, and to the next wavelength point 
/ -I- 1 . The A, B, C terms refer to the spatial points / - + 1, 



and rewrite Eq. (|48]l as 



(50) 



(51) 



-/new - A*S new + (A - A*)S old- 

This relation can be written as ( |Hamann|1987| l 

[1 _ A-(i _ e)] = 4 _ A*(l - 6)/oid, (52) 



respectively. Equations ( 25 i-(|44|i are nearly identical to those of 
Knop et al. ( 2009a| l except for the explicit Doppler factor f{s), 
which arises because the photon direction is measured by a dis- 
tant stationary observer here rather than by a comoving observer, 
as was the case in |Knop et al. ( 2009a I. We also clarifie d a prob- 



lem with C^:'' that was confusing in Knop et al. 



where 7fs = A5oid- Equation 52 is solved to obtain the new val- 
ues of /, which are then used to compute the new source function 
for the next iteration cycle. 

Mathematically, the ALI method belongs to the same family 
of iterative methods as the Jacobi or the Gauss-Seidel methods 



(2009a I. The for- (Golub & Van Loan 1989 1. These methods have the general form 



mal solution matrix is therefore identical m form to that shown 



in Figure 1 of Knop et al. ( 2009a i. 



M/+' ^Nx^ +b 



(53) 



An element of the source function vector Al is given by 
( |Knop et al.|2009al l 



(45) 



Note that all Doppler factors are explicitly handled by including 
them in the opacity as xi * /(■?)- 

From Eq. (35 1 the solution for the specific intensity for a 
given spatial point and wavelength reads 

rk _ /i Ddiag,':"\"' f >. rk , jjsub,*: r<: r/Saper.k .k 



.suh.kTk 
, y--isub,/; jk . 



.img,kjk , super,/: ,i 



^super,*^^; 



\,l+l> 



(46) 



for the iterative solution of a linear system Ax - b where the 
system matrix A is split according to A = M - N. For the 
ALI method we have M - l-A*(l-e) and, accordingly, 
N = (A - A*)(l - e) for the system matrix A - I - A(l - e). 
The convergence of the iterations depends on the spectral radius, 
p(G), of the iteration matrix G - M^^N. For convergence the 
condition p(G) < 1 must be fulfilled, this puts a restriction on 
the choice of A*. In general, the iterations will converge faster 
for a smaller spectral radius. To achieve a significant improve- 
ment compared to the A-iteration, the operator A* is constructed 
so that the eigenvalues of the iteration matrix G are much less 
than unity, resulting in swift convergence. Using parts of the ex- 
act A matrix (e.g., its diagonal or a tri-diagonal form) will opti- 
mally reduce the eigenvalues of the G. The calculation and the 
structure of A* should be simple to make the construction of the 
linear system in Eq. ( 52 1 fast. For example, the choice A* = A is 



2.2.1 . The operator splitting metliod 



Now that we have the formal solution, the fuU scattering problem 
can be solved by using operator splitting. The mean intensity 
is obtained from the source function 5^ by a formal solution 



best in view of the convergence rate (it is equivalent to a direct 
solution by matrix inversion) but the explicit construction of A 
is more time-consuming than the construction of a simpler A*. 

In the following discussion we use the notation of Hauschildt 
( [1992) 1 and [Hauschildt & Baron] ( [2006| ). The basic framework 
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and the methods used for the formal solution and the solution 
of the scattering problem via operator splitting are discussed in 
detail in Hauschildt & Baron (2006 ) and will therefore not be 
repeated here. We have extended the framework to solve line 
transfer problems with a background continuum. The basic ap- 
proach is similar to that of Hauschildt ( 1993[ ). In the simple case 
of a two-level atom with background continuum that we con- 
sider here as a test case, we use a wavelength grid that covers 
the profile of the line including the surrounding continuum. We 
then use the wavelength-dependent mean intensities 7^ and ap- 
proximate A operators A* to compute the profile-integrated line 
mean intensities / and A* via 



(p{A)JA dA 



and 



A* = J cl>(A)A* dA. 



J and A* are then used to compute an updated value for J and 
the line source function 

S = (l-£)J+eB, 

where e is the line thermalization parameter (0 for a purely ab- 
sorptive line, 1 for a purely scattering line). B is the Planck func- 
tion, B \, profile averaged over the line 

B = J (f>(A)BAdA 

via the standard iteration method 



[1 - A*(l - e)] /new = /ft - A*(l - e)/ 



old, 



(54) 



where /fs = A5oid- This equation is solved directly to obtain the 
new values of /, which are then used to compute the new source 
function for the next iteration cycle. 

We construct the line A directly from the wavelength- 
dependent A* generated by the solution of the continuum trans- 
fer problems. 



Given the form of Eq. 46 for the formal solution, the con- 
struction of the A* -operator can proceed exactly as described 
in [Baron & Hauschildt| ( |2004| . However, to conserve memory, 
we have implemented the A*-operator, retaining the spatial off- 
diagonal terms, but neglecting the off-diagonal terms in wave- 
length. We still find good convergence at considerable savings 
in memory (see below). 



3. Test calculations 

In this section we present the results of test calculations we per- 
formed to test the new algorithm in terms of accuracy by regres- 
sion testing. We compare these to results of the homologous case 
(I Baron et al.||2009b|l 



and to the spherical nonmonotonic cases 



( |Baron & Hauschildt|2004| . 

The test calculations were performed on Opteron CPUs run- 
ning Linux (Franklin and Hopper at NERSC), on Intel CPUs 
(Carver at NERSC, ICE2 at HLRN, and our own local Xserve 
clusters), and on IBM CPUs (PWR-4 and PWR-5). The code 
was compiled with Gfortran/gcc/g++, ifort/icc/icpc (versions 1 1 
and 12), and xlf/xlc/xlC and with NAG f95/gcc/g++. Using the 
varied compiler suites and CPUs allowed us to find numerous 
errors. 



3.1. Regression w/fh monotonic case 

Figure [T] shows the profile of the mean intensity J for homolo- 
gous flow and spherical symmetry. The test problem is similar to 
that of |Baron et al.| ( ,2009b) . These are the basic model parame- 
ters: 

1. An inner radius Ri„ - lO" cm and an outer radius /?out = 
1.01x10'^ cm. 

2. A minimum optical depth in the continuum t™" = 10""* and 
a maximum optical depth in the continuum t™" = 10'*. 

3. A gray temperature structure with Tinodei - 10^ K. That is the 
temperature solution of the spherical gray atmosphere prob- 
lem with effective temperature Tmodei (Mihalas 1978). 

4. An outer boundary condition = and an inner boundary 
condition 7^ = B^ for all wavelengths. 

5. For the initial wavelength the boundary condition is taken 
from that given by the 3D homologous calculation for ho- 
mologous tests and set equal to the Planck function B for 
non-homologous tests. 

6. A continuum extinction - C/r^, with the constant C fixed 
by the radius and optical depth grids. 

7. A parametrized coherent and isotropic continuum scattering 
given by 

Xc = £cKc + (I - ec)o-c (55) 

with < ff < I. Kc and cTc are the continuum absorption 
and scattering coefficients. In this work we have neglected 
scattering in the continuum. 

The line of the simple two-level model atom is parameter- 
ized by the ratio of the profile-averaged line opacity xi to the 
continuum opacity Xc and the line thermalization parameter e/. 
For the test cases presented below, we used = 1 and the line 
strength is given by F s xi/Xc - 10^ to simulate a strong line, 
with varying e; (see below). 

The test model is just an optically thick sphere put into the 
3D grid. The velocity at the outer radius was set to be relativistic, 
Vm-dx - 8x10"* km s"' . The calculation was performed on a spher- 
ical grid with 33^* spatial points, 33^ solid angle directions, and 
22 wavelength points. This and all calculations presented below 
were parallelized over characteristics and were run on parallel 
clusters. Ng ( |Ng| 1 974] |Auer| 1 987| l acceleration was used to sig- 
nificantly speed up the operator splitting iterations for cases with 
scattering. Figure|2]shows the line profile at the surface, here we 
compare the 3D monotonic calculation ( [Baron et aLpOOT) to 
the same calculation using the 3D arbitrary velocity algorithm. 
Because of the spatial resolution and the way that characteristics 
end at different points in a voxel (as opposed to always ending 
on a radial grid point in the ID case), it is better to compare 3D 
cases to 3D cases. 

Figure |3] shows the convergence rate of a monotonic flow 
case, with e - IQ~'^, with and without Ng accleration. Not only 
does Ng acceleration clearly improve the convergence rate, it 
gives almost the exact same errors as the same test setup using 
both the 3D and ID homologous algorithms. 

3.2. Sine velocity flow 

Here we again consider a spherically symmetric case with the 
same physical parameters as in § [3.1| but the velocity, while still 
radial is now given modulated by a damped sine wave. We have 
also introduced a line into the opacity with thermalization pa- 
rameter €i = 1. The thermalization parameter in the continuum 
is - 1. This case was calculated in ID in [Baron & Hauschildt[ 



5 



E. Baron et al.: Arbitrary Velocity Flows 




R (cm) 



Fig. 1. Comparison of monotonic flow case with spherical sym- 
metry solved using the full arbitrary velocity field method (red) 
to the well-tested ID solution (black). The comoving mean in- 
tensity is plotted at each point in the computational volume. The 
agreement is at the 1% level. 



( |200 4). Figure |4] shows the velocity as a function of radial opti- 
cal depth in the continuum. The spatial grid was 129 x 65 x 17 
and the solid angle resolution was set to 512^. The maximum ve- 
locity is V X 10, 000 km s"' and the total number of wavelength 
points is 226. Figure |5] shows the comoving mean intensity 
from each surface voxel as a function of wavelength A. The ID 
result is plotted (the green curve) and the agreement is good to 
the sub-percent level. The variation in the velocity leads to an 
asymmetric line profile even in the comoving frame. But since 
this test has such a high spatial resolution, it requires a signifi- 
cant amount of memory per process. We explored the effects of 
lower spatial resolution, 129 x 17 x 17, while keeping the high 
soUd angle resolution 512 x 512. Figure |6] shows the mean in- 
tensity as a function of wavelength for the outmost part of the 
sphere. For computational expediency we set the line thermal- 
ization parameter e/ = 1. The spread in the results of 1.5 % is 
indicative of the accuracy of these calculations, whereas the de- 
viation of the envelope of the solution is indicative of the low 
spatial resolution of this calculation. While modern nodes may 
have 24 CPUs, but only about 1 GB of RAM per CPU, it is 
unfortunate that the evelope of the low spatial resolution calcu- 
lation is offset from the ID result. This shows that not only solid 
angle resolution is important, but spatial resolution is quite im- 
portant as well. Figure |7] shows the same calculation with the 
spatial grid enhanced to 129 x 65 x 65 and the solid angle grid 
reduced to 32 x24 ( jHauschildt & Baron|2010| l. As can be seen in 
Figure |7] the reduced solid angle resolution increases the spread 
by roughly a factor of two, but now the 3D solution envelopes the 
ID nonmonotonic result. With this higher spatial resolution, the 
spread can be further reduced by increasing the solid angle res- 
olution, which shows near-perfect weak scaling, and thus does 
not increase the wallclock time required for these calculations 
(although it does require more CPUs). 



3.3. Example of radial nonhomologous flows 

The previous tests were all totally spherically symmetric, with 
radial variations in the velocity field. We now assume the flow to 
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Fig. 2. Comparison of monotonic flow case with spherical sym- 
metry solved using the full arbitrary velocity field method 
(green) to the 3D monotonic flow solution (red). The comov- 
ing mean intensity is plotted at the surface of the sphere. The 
test cases are identical, Vmax = 8 x 10** km s"', 129 x 17 x 17 
spatial voxels, 1 14 wavelength points, and 256 x 256 directions. 
The thermalization parameter in the line is e = 10"^ and the line 
strength is F = 100. The lower panel shows the relative differ- 
ence, 6Ja/Ja, as a function of wavelength. 



be radial and azimuthally symmetric, i.e. 

m = p(0)q(r)t,. 

For this case, we have 

h - J3 - pqn,; 

1 



and 



(n ■ V)y6 = pn,q'(r) + -p'{e)qng, 
r 



1. 



(56) 

(57) 
(58) 



n ■ V(n ■ ;8) = -[p'qrirng + p(q'r)nf + pqil - <)], (59) 



where we made use of 
r — n,., 

and 



rig . 

— , (p - , 

r r sm 



rir - ngO + sm 



(60) 



(61) 



Here an over dot implies ^, e.g., r = ^, etc., and a prime de- 
notes differentiation with respect Xo ji - cos 9. 
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Fig. 3. Comparison of the convergence rate of a monotonic flow 
case with spherical symmetry with the thermalization parameter 
e = 10""* with and without Ng acceleration. The Ng acceleration 
is identical to that produced by the pure homologous 3D module. 
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Fig. 5. Line profile for the damped sine-wave velocity case with 
high spatial - 129, ng - 65, and = 17, and angular resolu- 
tion «n - 512 X 512. The comoving mean intensity for each 
voxel on the surface (there are 65 x 17 = 1 105 of them) is plotted 
as a function of A (red lines). The green curve is the ID result 
using the methods of Baron & Hauschildt ( 2004| l and Knop et al. 
(|2009b! ai Here the scattering fraction is e = 10 and the line 
strength is T = 100. 



symmetric spherical functions) 



P(0) = 2^«^«^^^' 



(64) 



with Po(ij) = 1, PiQu) = iJ, and 'P2(m) = {O/j^' - 1), etc., where 
fi = cos 6. We consequently obtain 



p'(0) = ^^«^«^^) 



(65) 



Fig. 4. Velocity profile used for the nonmonotonic velocity test 
with a damped sine wave. 



Inserting Eqs. ( 57 i-(|59]l into Eq. (19i, we find an analytic 
expression for a(i): 



a(s) — 



p'qrirrig + p{q' r)n; + pq(l - nf) 



r(l - pqrir) 



-J^pq 



prirq'ir) + -p'{0)qnt 
r 



(62) 



where the spherical components of the unit vector fi = 
(sin 0„ cos 0„, sin 0„ sin <p„,cos 9n) are 



rir — sm6sm0„cos(<p - (p„) + cosOcosOn, 
rig - cos0sin0„cos(0 - 0„) - sinflcos^n, 
— cosd„. 



(63) 



For p(ff), we use an expansion in terms of Legendre polyno- 
mials (which form a complete and orthonormal basis for axially 



with !Pj - - sin 0, and '^2 ~ ~ 2 ^™ (i^o^e that the case 

c„ - C()5° degenerates to homologous flow). We can use this 
finite expansion in terms of the Legendre polynomial 'P„i0) to 
construct azimuthally symmetric jets. 

As a simple example of nonspherically symmetric flow, we 
run a test case where 



(66) 



i.e. we take q{r) - c^r (cq is a constant), which simplifies a(i) to 
be 



a(s) = Co 



P+ p HrUg 



1 - rrirp 
Furthermore, we assume 



j^Piprir + p'rig) 



p(0) = i + ^^2(/^) = i + M^cos^e-l 



(67) 



(68) 



i.e., we perturb the homologous flow by adding a second-order 
Legendre polynomial with perturbation coefficient 0.5. We show 
the resulting mean intensity plot J{6, (p) at the boundary R - /?out 
in Fig. [8] We obtained perfect azimuthal symmetry as expected 
from the form of in Eq. ( 66 1, although it was not explicitly 
enforced. Furthermore, we also recovered the symmetry with re- 
spect to the north/south pole (i.e., symmetry under a reflection 
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Fig. 6. Line profile for the damped sine-wave velocity case with 
lower spatial - 129, ng - \1, and - \1, and higher angular 
resolution nn = 512x512 than in Figure |5] The comoving mean 
intensity Jx for each voxel on the surface (there are 65 x 17 = 
1 105 of them) is plotted as a function of A (red lines). The green 



curve is the ID result using the methods of Baron & Hauschildt 
( |2004| l and |Knop et alT] (|2009b'a'). Here the scattering fraction 
is e = 1 for computational expediency, and the line strength is 
r = 100. With the very low spatial resolution, the 3D result 
no longer envelopes the ID result (which by assuming spatial 
spherical symmetry has essentially infinite resolution in ng, and 



Fig. 7. Line profile for the damped sine-wave velocity case with 
higher spatial resolution, but lower angular resolution than in 
Figure [5] The spatial resolution is - 129, ng - 65, and n^j = 
65, andthe angular resolution is reduced to na = 32 x 24. The 
comoving mean intensity 7 ) for each voxel on the surface (there 
are 65 x 17 = 1105 of them) is plotted as a function of A (red 
lines). The green curve is the ID result using the methods of 
Baron & Hauschildt] ( [2004l l and |Knop et~l. ( 2009b a ). Here the 
scattering fraction is e = 1 for computational expediency, and 
the line strength is F = 100. 




with respect to the equatorial plane, ^ n — 0), which is char- 
acteristic for Legendre polynomials of even order, see Eq. (j68|. 
Figure |9] shows a longitudinal slice of the comoving intensity, 
which shows that the comoving mean intensity varies by roughly 
a factor of ten from the pole to the equator. This is just due to the 
effect that the continuum level varies with the maximum velocity 
(Baron et al. 2009b| ). 



3.4. Checkpointing 

We implemented a checkpointing scheme that allows for restarts 
and new starts with higher resolution in solid angle. We simply 
write out (using stream I/O) the A* operator once it has been 
calculated, and the value of J at each voxel after each ALO it- 
eration. Since these calculations require significant numbers of 
processors, which may go down, or calculations may run out of 
time, this allows us to perform perfect restarts. We checked that 
restarting from the checkpoint files works perfectly and that the 
calculations are converged in a single iteration when restarting 
from a converged checkpoint file. We write A* and J to sepa- 
rate files, since that allows us to read just J, and then recalcu- 
late A* if we restart with a higher resolution. We checked this, 
and it allows us to use a low-resolution calculation (performed 
on fewer processors, for example) and then converge the higher 
resolution calculation with roughly half as many iterations for 
each resolution increase. Table [T] shows the rate of convergence 
for a test with homologous velocity fields (chosen for compu- 
tational expedience), v,nax = 10^ km s ' and e = 10"^. While 
the restart result still requires several iterations, this scheme al- 
lows for some speedup. The small e and high maximum velocity 
make this test particularly demanding. The calculations were all 
restarted such that Ng acceleration cannot begin until the fourth 




Jmin Jmax 

Fig. 8. Comoving mean intensity J{6,(p) plotted at the surface 
R - Rout for case where jS = CQrp{6)tr, and p{6) - \ + ^Pii^i) = 
1-1-5(1 cos^ 6- 5). The velocity flow is radial, but not spherically 
symmetric to approximate a jet-like flow in the +z direction. 



iteration, and thus the restart calculations converge more slowly 
than they could. The scheme could indeed possibly be made 
more efficient by keeping the previous three values of /, so that 
the restart iteration could immediately use Ng accleration. 

4. Conclusion 

We presented algorithm strategies and details for the solution of 
the radiative transfer problem in 3D atmospheres with arbitrary 
wavelength couplings. 
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10 20 30 40 50 60 70 80 90 

Fig. 9. Comoving mean intensity J{6,(p()) plotted at the surface 
R - 7?out for a fixed (arbitrarily chosen) value of <p - (pQ. The 
variation of nearly an order of magnitude from the pole to the 
equator is due to the fact the the continuum level is a function of 
the velocity. 

Table 1. Convergence rate for homologous test with Vmax - 
10^ km s"' and e - 10"^, starting from the previous test. The 
first test with 8x8 solid angles is converged from the beginning. 
The spatial resolution is held fixed and A* is recalculated at the 
beginning of each new resolution run. 

n Number of iterations 

~8xl 29 

64x 64 21 

128x 128 15 

256x256 13 



Future possible applications are the velocity profiles of cool 
stellar winds, treatming partial redistribution, calculating radia- 
tive transfer in shock fronts like in accretion shocks, calculating 
quasar jets, and general relativistic situations such as a rotating 
black hole. 
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